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C.I 



ABSTRACT 



This thesis models the effects in a fluid medium of a plane wave that has impinged upon a 
reinforced plate. The wave equation for pressure, and the equation of a thin plate combined with 
other equations are coupled at the interface between the fluid and the thin plate. The actual modeling 
is accomplished in a fortran computer program written to run on the Naval Postgraduate School’s 
main frame computer. The program uses extensive finite differencing on a domain, assumed to be 
a small section of an infinite interface between the fluid and the plate, to simulate the deflection of 
the thin plate and the pressure disturbances in the fluid medium. To accomplish this, each of the 
above equations will be scaled or nondimensionalized. Additional finite differencing is explained 
which covers the special cases for the side boundaries of the fluid domain, and the artificial boundary 
created to model infinity. Different beam spacing is explored for its effect on the magnitude of the 
propagating modes of the scattered pressure wave. 
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I. INTRODUCTION 



This thesis models numerically, the effects of a pressure disturbance incident on 
a fluid loaded reinforced plate using finite difference methods. The program produced 
explores the reaction of a plane wave striking the plate (hull of a ship) at different 
frequencies, for different plate lengths and for various spacings of the reinforcing beams 
(or members). 

Figure 1 shows the problem being modeled. The figure shows a cross section of 
a reinforced flat plate with water on the top side of it and a vacuum below it. For this 
problem, it is assumed that the plate is infinite in the x direction along the horizontal 
axis. Since the plate structure is periodic in space, the model will only look at one 
section of this infinite plate. The behavior at other sections is assumed to be the same. 
We will use three fundamental equations in the model. The first equation is the wave 
equation, the second is the linear-inviscid-force equation, and the third is the equation 
of motion for a thin plate. The coupling of these equations at the fluid-plate interface 
makes the problem very difficult to solve analytically. 

Chapter II derives the relationships among these three equations as they are used 
in the model. Assumptions used in the program are spelled out in this chapter. Chapter 
II also includes the steps necessary to scale each equation to facilitate its use in the 
computer program. 
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Chapter III explains the operation of the FORTRAN computer program used to 
model the problem. The program takes as an input the incident and specularly-reflected 
sound waves and determines the scattered pressure due to the fluid-structure coupling. 
To make best use of the Naval Postgraduate School’s mainframe computer the main 

i 

program is divided into a number of smaller subroutines. These are also described in 
Chaptdr III. 

Chapter IV discusses the results of the program runs. It also reviews the 
weaknesses in the model and model limitations. 
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Figure 1. Cross Section of Hull 
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II. THEORY 



A. THE WAVE EQUATION 

It is known that sound travels through seawater, or any fluid medium, as a pressure 
wave. This is possible because seawater is a compressible fluid which is affected by 
pressure disturbances. To approximate the movement of sound in a fluid medium we 
have used the linearized, lossless wave equation for pressure as derived in Kinsler, Frey, 
Coppens and Sanders [Ref. l:p. 104], 






1 




( 2 . 1 ) 



where Cf is the fluid sound speed and p^ is the total pressure. 

The pressure wave can be split up into various parts, each of which is a solution 
of the linear wave equation. In this case we want to be able to solve for pressures that 
would be present in the fluid which modify the pressure resulting from specular reflection 
of an incident plane wave from an acoustically hard surface. This known specularly 
reflected pressure wave is combined with the computed scattered pressure wave to make 
up the total reflected pressure wave. 

Symbolically, we have; 
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p^=p'^p’^+p^, 



( 2 . 2 ) 



where is the total pressure, p' is the incident pressure, p'^ is the specularly reflected 
pressure, and p® is the scattered pressure. The total radiated pressure is a superposition 
of the reflected and scattered pressures. 

To eliminate the specularly reflected wave from the equations of motion we apply 
the boundary condition 



waves must cancel at the surface of the plate. The vertical coordinate Y is equal to zero 
at the surface. 

We assume that the source of the incident pressure is far from the surface so that 
a plane wave approximation is adequate to model the incident pressure, p‘. 

The time harmonic incident and specularly reflected plane waves are of the form: 




^ ^ y-O 



(2.3) 



In this case the normal derivatives of the incident and specularly reflected 




(2.4) 



4 



and 



P 



\x,y,t) = e 



zu( -0 



(2.5) 



where 






COS0^ 

-sm9, 



( 2 . 6 ) 



and 






cosQj 

sin0^ 



(2.7) 



are the unit direction vectors of the two plane waves. 

Both the incident wave, Equation 2.4 and the reflected wave, Equation 2.5 must 
satisfy Equation 2.1. Therefore, because of superposition, p* must also satisfy the wave 
equation. 






1 

cf dt^ 



( 2 . 8 ) 



B. THE LINEAR INVISCID FORCE EQUATION 

The linear inviscid force equation, taken from Kinsler, Frey, Coppens and Sanders 
[Ref. l:p. 104], is basically the inviscid Navier Stokes equation in which nonlinear terms 
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have been neglected. It is valid for acoustic disturbances of small amplitude, and can be 
written 

^ Ot 

where Pf is the constant equilibrium density of the fluid, and u is the velocity of the 
fluid. 

Since we assume no cavitation of the fluid or penetration of the fluid at the fluid- 
plate interface, the component of the fluid velocity normal to the interface and the plate 
transverse velocity are equal at the fluid-plate interface at all times t. Letting w represent 
the transverse displacement of the plate (y direction), we have 

„ ,, = 0 . ( 2 . 10 ) 

‘ Sy 

Through use of Equations 2.2 and 2.3 the total pressure p^ can be replaced by the 
scattered pressure p®: 

= a, y^o, ( 2 . 11 ) 

C. EQUATION OF MOTION FOR A THIN PLATE 

To properly understand what is happening on and within the thin plate, a reference 
system must be assigned. In Figure 2 we see that the vertical plane will be designated 
the y-axis in which the thin plate experiences transverse displacement w. The horizontal 
plane will be the x-axis which has longitudinal displacement u. To set up our equations, 
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we assume that the centerline (axis) of the beam stays perpendicular to its cross section 
during bending. 
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Figure 2 Plate Reference System 



As put forth by Kinsler, Frey, Coppens and Sanders [Ref. l:p. 58], Hooke’s law 
states that if the strain is small the stress is proportional to it. The stress on the plate is 
described by the equation 



o = -y 



E 



( 2 . 12 ) 



where: 

<r„ is the stress on the plate, E is the Young’s modulus, and u is the Poisson’s ratio. 
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Consider a small section of the plate of area AxAz with moments and forces applied 



as shown in Figure 3: 
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Figure 3. Forces and Moments 

Note that there are no forces or moments applied which have z dependence. The 
motion therefore is assumed to be two-dimensional. Since the plate is not rotating, a 
force balance is performed on the moments and moment arm, neglecting terms of order 
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Ax^, giving us 



dM 

dx 



(2.13) 



where: 

M = bending moment on the element, 

q = load on the plate element (uniform for a normally incident wave), 
and F = flexural force applied to the element. 

Equating forces in the y-direction applied to the section of plate yields 



where: 

S = cross sectional area of the plate, 
p, = density of beam section, 
h = plate thickness. 

Dividing out the AxAz terms yields 




(2.14) 




(2.15) 



Using the definition of the moment 



h 




(2.16) 
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where the moment of inertia is found from 
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we obtain the equation for deflection of a plate 




(2.18) 



The loading of the thin plate is due to pressures exerted by a fluid in contact with 
the thin plate. Rewriting Equation 2.18 including the pressure terms yields: 



is the flexural rigidity and h is the plate thickness, Junger and Feit [Ref. 2:p. 194]. 

D. THE SCALED EQUATIONS 
1. The Scaling Factors 

In order to use the three governing Equations 2.1, 2.9 and 2.19 in the finite 
difference code, they must be scaled or nondimensionalized. The steps are as follows. 
The independent variables x and t are scaled through use of the fundamental length L 
and the wave frequency w, respectively. 




(2.19) 



in which 




( 2 . 20 ) 



and y = — , 

L 



( 2 . 21 ) 



x = 



I’ 



where x and y are the scaled spatial coordinates, and L is the one half the distance 
between stiffeners measured in the x-direction. 

The equation for the scaled time is 



t = or, 

where « is the frequency of excitation of the incident plane wave. 
We then have 



( 2 . 22 ) 



a _ a _ 1 a 

dx d(Lx) L dx ’ 



( 2 . 23 ) 



and 



a _ a _ _a 

ar 7 ^ dt 

d(-) 

CO 

Additionally, the dependent variables w and p must be scaled as follows. For 
the displacement of the plate (in the y-direction) we put 



- w 




( 2 . 25 ) 



where w is the scaled vertical displacement of the plate and W is the displacement 
scaling factor. For pressure, put 
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(2.26) 




where p is the scaled pressure and P<, = Lw^pfW is the pressure scaling factor. 

The wave numbers for waves propagating in the fluid, and along the plate are 
given as follows, 




(2.27) 



and 



(2.28) 

respectively. 

Additional parameters resulting from the scaling analysis are 



I’ D 



kj 



CO 



(2.29) 



kl 



where Q is the squared ratio of fluid to plate wave numbers and 5 )q, the coincidence 
frequency of the plate, is defined by 



‘^ 0 = 1 - 



Jicf T 



/ l2 



D 



We also define the dimensionless quantity 



(2.30) 
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6 = - 



P// 



(2.31) 



where e is a dimensionless quality which is a measure of the fluid loading on the plate. 
2. The Wave Equation 

The wave equation governing the scattered fluid pressure was found to be: 



cJ dt^ 



(2.32) 



Scaling the above equation yields: 



1 1 

Cf L L 



(2.33) 



which reduces to the desired non-dimensional form, 



- - . 21 . 2 - 



(2.34) 



3. The Linear Inviscid Force Equation 
The linear inviscid force equation 



P/ 



_ dp ‘ 



dt^ dy 



in scaled form reduces to 



(2.36) 



W--= -D-. 
1 1 



(2.37) 
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4. Thin Plate Equation 



The dimensional equation of flexural rigidity is: 




(2.38) 



By substituting in the scaling factors it becomes: 



D 



(2.39) 



or, after more simplification, it becomes 




(2.40) 



This is the scaled form of the thin plate equation that will be used in the 



computer program. Note that at y = 0, p*^ equals p'. Hence, p’’’ has been replaced with 



5. Stiffeners 

In the present analysis the stiffeners will be treated in an idealized way. The 
end conditions we will assume are imposed at the point where the stiffeners are located 
will be equivalent to those of a clamped plate for which the boundary conditions are 
w = 0 at X = ±1, and w^ = 0 x = ±1. 

E. THE FLUID DOMAIN 

As stated before, it is assumed the plate on which the sound wave impinges is 
infinite. That assumption makes it possible to establish the following conditions. The 
boundary conditions on the section of plate we are looking at will be periodic. That is 



P + 2p‘. 
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to say, for normal incidence of p', displacement on the plate at x = -1 will be the same 
as at X = 1. For arbitrary angles of incidence, the pressures at x = ±1 differ only by 
a scaling factor determined by the angle of incidence. Figure 4 and the following 
equations will show how the time harmonic solution is dealt with in the problem. 
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Figure 4. Side View of Fluid Domain 



To understand how the fluid pressure is disturbed, the wave equation is considered 
for which a prescribed pressure is applied from the fluid to the surface of the plate. This 
uncoupled problem can be solved analytically by the method of separation of variables. 
The scaled wave equation is 
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(2.41) 



P +P =a^P . 

XX yy tt 



Applying periodic boundary conditions results in 



p(-ij)=;>(i,y) and ;>;,(-i,50=;>;,(iJ)- 



(2.42) 



The forced pressure at the plate surface is given by the boundary conditions 



p{xfi,t)^g(x)e''^ where g(-l)=g(l) and g^(-l)=g,(l), 

where: 

g(x) = prescribed pressure applied to the fluid by the plate. 

Assuming a time harmonic solution with equivalent frequency dependence to that 
of the forcing function, the scattered pressure can be expressed as 

By substitution of Equation 2.44 into Equation 2.41 we obtain the reduced 
Hehmholtz equation for 4>(\,y), 



(|)i7+%+a'<t>=0, (2-45) 

with boundary conditions: 

4>(-l J) = 4>(1,5’) and <t>^(-l,y)=4>;t(l,5') (2-46) 

and 
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(2.47) 



Henceforth we will be dropping the bars from all equations, giving us the 
following, 



4>„+4’^.y+«^<}>=0, 


(2.48) 


4>(-l,y)=(})(l,>-), 


(2.49) 


4>.(-i>y)=4>.(ij). 


(2.50) 


<l>A-hy)=Hiyy 


(2.51) 



To solve for <f> the standard method of "Separation of Variables" is applied to the 
equation. First we write (/> as a product of functions of one variable each, 

4>=A:(A:)y{». (2.52) 

Upon substitution into the Helmholtz equation and separation of variables, we find a 
solution 

X^(x)=C/‘”’'*. (2.53) 

When solved for Y, two answers result: 

when a>/m, (2-54) 



and 
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when cc < AJ 71 . 



(2.55) 



These solutions are chosen to avoid the physically implausible solution that the scattered 
pressure as y approaches infinity will blow up represent an incoming wave. 

Combining our solutions we have 

- N 

4>=52<})^(x,y)=5]) ^ ( 2 . 56 ) 

n=0 n =0 n=N*l 

Note that as y approaches infinity the decaying terms become exponentially small. Now 
from our boundary condition at y = 0, a„ can be found such that 

“ 1 * 

g(x)=^ with a^=— fg(x)e'‘"”^clx. (2-57) 

It is important to note that in the numerical solution of the fully coupled problem, 
the radiation waves must be properly modeled as y approaches infinity. The numerical 
domain cannot be semi-infinite in extent in the y direction. It must be truncated at some 
value of y, which we call y,,, where a boundary condition must be imposed to ensure 
proper modeling of the radiated pressure. 

F. THE PERIODICITY PROBLEM 

To ensure that a wave with any angle of incidence can be correctly modeled by the 
computer program the periodicity must be worked out. If we were to have a wave 
normally incident on the plate we would expect a solution for the pressure which is 
periodic in x to be 
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(2.58) 



4>=E 

rt =0 

To account for an arbitrary angle of incidence, Equation 2.58 needs to be modified 
in the following manner. A correction factor which takes into account the phase change 
at X = -1 vice x = 1 must be added. The modified equation for the scattered pressure 
is, 



rt=0 

where is the angle of propagation of the incident pressure wave. 
Continuing with this argument one can write 

Y =-^COS0,+rtTt. 

* n i 

Substitution of this x dependence into the Helmholtz equation yields 




dy^ 






Substituting these relationships back into our equation for </> gives us 

N -M-\ 

*=E E *E 

n-~M -00 //+1 

where 



(2.59) 



(2.60) 



(2.61) 



(2.62) 
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(2.63) 



\/ y !^ k<y„ 



and in which the first sum includes all values of n for which k > 7„ and the latter two 
sums account for the remaining modes. 

G. FINITE-DIFFERENCE MODELING 
1. The Taylor Series 

Refer to Numerical Solution of Partial Differential Equations: Finite 
Difference Methods . Chapter I, by G. D. Smith [Ref. 3] for a more complete explanation 
of finite difference methods. 

All three governing equations will be center differenced with truncation errors 
of order (Ax^, At^). 

Finite difference approximations for u„ u„, u„„ are 




(2.64) 




+0(/i2) 



(2.65) 



and 




(2.65) 



where Uj = u(Xj). 



20 



2. Differencing of the Wave Equation 

We have already shown that the wave equation is of the following form 



P.. -Pyy = 



(2.67) 



and that we have the following conditions, 



a:. = jAx, yj = j^y, = nAt, 



( 2 . 68 ) 



where i, j, and n are integers. This allows us to describe the pressure as follows, 



Pij = Ax = Ay = h. 

So by substituting into Equation 2.66 we arrive at the final form. 



(2.69) 



P 



n*l 



— — IP", +P", •-4P"+P" ,+P" A+2P" -P”'\ 

^ 1-1, y ij ^ ij*\ ^ ij ^ ij * 






(2.70) 



3. Finite Differencing of the Linear Inviscid Force Equation 

We have already shown that the linear inviscid force equation is of the 
following form 



dp 

-w, =-^. 

" ay 



(2.71) 



This equation is applied at y = 0 or when j = 0 in the numerical domain. 
The displacement w does not depend upon y, so we introduce W" = W(Xj, t„ ) and write 
the differenced form of the equation as: 
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(2.72) 



' 2Ay 

This equation is solved for the P"; .1 which is exterior to the numerical domain. 

It can then be substituted into the differenced form of the wave equation to solve for 
P“'^*io > the pressure on the plate fluid interface at time level n + 1 solving for Pj., we 
have 



+ (2.73) 

Ar^ 

Note that is first needed to solve for P"j .| which means a prescribed order must be 
followed at the fluid plate interface. First W"^‘j is found, the P"i,.i, and finally P"‘^'i.o- 

4. Finite Differencing the Plate Equation 

The last of the major equations that we must put in the form of finite 
differences is the thin plate equation. The scaled thin plate equation is: 



W, 



XXXX p 



-kyWj; ■ 



(2.74) 



Differencing this equation results in an expression for the displacement W at 
time level n + 1 in terms of the preceding two time levels and the incident plane scattered 
pressure at time level n. 
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5. Finite Differencing Side Boundaries 

When differencing the wave equation at the side boundaries for the pressure 
we have a similar problem to that found in applying the wave equation at the fluid plate 
interface. The problem is that there are values of the pressure needed which are exterior 
to the numerical domain. Periodicity will be used as an additional condition to eliminate 
this problem. We have the scattered pressure in the form: 

<2.76) 

«=0 

What we are doing here is taking the pressure at the node i = M + 1, x = Ax 
and replacing the pressure i = -M + 1 and x = -1 +Ax. From Equation 2.76 we obtain 

<2-28) 

and we also have: 

( 2 - 79 ) 



Now we can make the following two statements: 



(2.80) 








(2.81) 



Combining the previous statements we get: 




-2iJtcos0 




(2.82) 



An equivalent argument used at the left boundary yields: 




likcosO 




(2.83) 



6. Finite Differencing the Radiated Boundary. 

The final edge of the numerical domain is the artificial boundary which must 
allow free passage of radiated pressure disturbances. It must be far enough away so that 
the evanescent modes can be neglected. To see how the boundary condition is derived 
we start with the scattered pressure p* and try to derive an analytical expression for the 
pressures normal direction at the boundary, as is done in Scandrett and Kriegsmann 
[Ref. 4]. This type of radiation boundary condition is often referred to as a non-local 
radiation boundary condition. 



From previously we have 



N 

p^(x,y,t) =e evanesent modes). 



(2.84) 



where 



24 



= -ikcosdj+nn, 



( 2 . 85 ) 



and 

V^p^=k^p,^. 

Additionally in our problem we have 

k=k^L and k^^—, 

but as y approaches infinity p® is equivalent to only the radiated modes, 
time dependence into the a„ coefficient such that, 

P"- E a.COe'"-”'-'’, 

n = -M 

and in the following steps use orthogonality to find an(t)e'^"^. We have 

^ n ^ 

j e ^(x,y,t)dx = ^ '^'^e'^'^dx 



lr[-iJ{xos0,+rt7i -{-ikcosdf’^mn)] 

=e ' ' 

^ LX[nTt-niTt] 

We can therefore write 

where | = dummy variable of integration. 



( 2 . 86 ) 



( 2 . 87 ) 



( 2 . 88 ) 

We include the 



( 2 . 89 ) 



( 2 . 90 ) 



( 2 . 91 ) 
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(2.92) 



a„(t)e = ^fe 
^-1 



Noting the similarity to zero in the following equation 



N 



y [i..p 

.Wdy ’ 



Which can be rewritten as 



(2.93) 



0. 

dy 2 „..m ot 



(2.94) 



Remember, these equations pertain to the artificial boundary which is modeling infinity. 
For this problem y is equal to ye at the artificial boundary. 

When we start to finite difference Equation 2.94, the artificial boundary and 
use the trapezoidal rule of integration we obtain; 



pm pm . N IX 

E P> E -Pw'))- -0 <2-«) 

2n 4At„..Af /.-/X 



in which we take: 



s,= 



- l = ±IX 

2 

1 l*±IX. 



(2.96) 



Rewriting Equation 2.95 we have 
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(2.97) 



P” -P 
2h 



IX 



N 



■ " E M E 



4A<, 



K'’w'-fw‘)-0 



« = -M 



or 



-P 



iJY-\ 



IX 



2h 



h-IX 



(2.98) 



Where the square matrix A is defined as 



A 






2h 

4Ar 



N 



6,E 






Combining the above equations gives us; 



(2.99) 



IX 



PCn.^-PCiy-r^Y: A,iP:,;^-p-^)-0 



^m-lx 



(2.100) 



/=-/x 

From the differenced form of the wave equation applied at j = lY (the artificial 
boundary) we have: 



pB + l _ _p>'”l 

^i,IY ~ ^i,IY 



All 



[p-. 



+p" +p" 






4Af^ 




( 2 . 101 ) 



Combining Equations 2.100 and 2.101 and eliminating P\iy+i we obtain: 
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( 2 . 102 ) 



= -p" rp" +p" +2P” 1 

/=-/x 



k}h^ 



Now we need to put this equation into a form which makes it possible for us 
to see the matrixes needed to solve the problem. First we write; 



and 



A ^ /=-/x 

A 2 ^ A 2 

E -4,/w 

/i /: /--/X /i k 



(2.103) 



We define 



6 =fl 

lO if i*l 



(2.104) 



A,,-6,+ A., 



(2.105) 






B , = -6;,+ A., 



(2.106) 



Equation 2.102 can then be written in matrix form as 



Ap;;.'* = Bp;;:'+c 



(2.107) 



where 
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(2.108) 



h k h k 



Then the equation for the pressure vector at the next time level is given by: 



P^"'=A''[BP”'’+q (2.109) 

The pressure at P"^‘i,iY is then assigned by taking the i*^ component of the 
vector In solving the matrix Equation 2.107 we use the UNPACK library 

Gaussian elimination programs (CGESO and CGESL). 
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III. DESCRIPTION OF PROGRAM 



A. GENERAL 

The computer program written for this thesis simulates the effects of a pressure 
disturbance on a fluid loaded reinforced thin plate. A detailed description of the theory 
behind the program can be found in Chapter 11. A pressure wave in a fluid medium is 
moving towards a reinforced plate wliich is flat, and is supported by stiffeners which are 
at a constant spacing on the non-fluid side of the plate. The non-fluid side is assumed 
to be a vacuum. When the pressure sound wave hits the plate there will be some (a very 
small amount) distortion of the plate in the vertical direction. This displacement will 
cause a disturbance in the fluid medium. The scalded pressure amplitude for different 
frequencies and beam spacing will then be plotted and the magnitude of the propagating 
modes will be recorded. 

The program stores the time dependent wave equation with time harmonic forcing 
by allowing transient effects to radiate away from the numerical domain. When a steady; 
state has ben achieved the code automatically stops the calculations. 

B. SUBROUTINES 

To make effective use of the Postgraduate School’s mainframe computer, the 
program has been broken down into a number of smaller programs (i.e. subroutines) 
which can be more efficiently programmed and run. 
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The main program is used to coordinate the subroutines and to pass information 
from one subroutine to another. The following are descriptions of the subroutines. 

1. Subroutine INIT 

The INIT subroutine initializes the computer program. In this subroutine we 
set up our two main matrices. The first of these is the "P" matrix, which represents the 
fluid pressure in the fluid domain. The second is the "W" matrix, this represents the 
displacement of the plate surface which bounds the fluid medium. 

Although we have the ability to change any number of parameters in the 
program, we will concentrate on two. These are the beam spacing, and the frequency 
of the pressure wave impinging on the plate. Both of these values are assigned in this 
subroutine. The INIT subroutine sets all of our counters and subscripts to zero at this 
time and initializes the matrices which will be needed to model the artificial boundary. 

2. Subroutine DOMAIN 

The DOMAIN subroutine calculates the pressure at all nodes in the interior 
of the fluid medium at a subsequent time level in terms of the preceding pressure values. 
The subroutine calculates the pressure at the next time step by use of the differenced 
wave Equation 2.70. 

The equation above must be modified to model what is happening at the side 
boundaries of the fluid domain. As was stated earlier, our model is looking at a small 
section of an "infinite" fluid plate interface. This allows us to use the pressure in the 
next section to model what is happening in the section in which we are interested. This 
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is necessary since otherwise we would have to use "undefined" points in our differenced 
wave equation at the edges of the fluid domain. 

3. Subroutine FLPL 

The FLPL subroutine models the fluid-plate interface along the top of the 
reinforced thin plate. Again, with finite differencing and the assumption that we are only 
looking at a small section of an infinite interface, we are able to model the effects of the 
pressure wave along the entire fluid-plate interface. This subroutine is also where we 
introduce our driving function into the computer program. This function is what starts 
things moving towards a steady state. Although this is not the most complicated 
subroutine in the entire program it is the one which has proved to be the hardest to make 
work as predicted. 

4. Subroutine ARTF 

The ARTF subroutine is the most complicated of any subroutine in the 
computer program. It models the artificial boundary of the fluid medium at infinity. 
Again, this is accomplished with finite differencing and matrix operations. 

5. Subroutine ALPHA 

The ALPHA subroutine calculates the magnitude of the radiating modes of the 
pressure wave. This data is gathered in the results section of this study. 

6. Subroutines SHIP and FINN 

The SHIF subroutine moves the entire problem to the next time level. This 
is done after the other subroutines have completed all the required calculations needed 
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to satisfy our domain and all of our boundary conditions. After exiting the SHIP 
subroutine the entire program is run again with the exception of the INIT subroutine. 

The FINN subroutine is a check to make sure that we have reached a steady state 
situation before we take any measurements of the plate displacement or of any fluid 
disturbances. Figure 5 is a graphical representation of the displacement convergence 
factor approaching zero. This was arrived at with a assigned frequency of 3500 Hz and 
L = 1.0. 

7. Subroutine CLOSE 

The CLOSE subroutine is used to send the data generated in other parts of the 
main program to the plotting devices. 

C. PROGRAM VERIFICATION 

Verification of the program was accomplished by the following method. We solved 
the waveguide problem (section 2E) for a known value of g(x). We let g(x) = sin ttx and 
e""’'. This was accomplished by replacing the FLPT subroutine with a subroutine named 
TEST. As before the program was run until a steady state solution was reached. These 
results were checked against analytical solutions and found to compare very well. For 
g(x) = sin TTX we graphically reproduced the expected results. For g(x) = e'” we 
calculated the coefficients and found Oj ~ 1 and all other a„’s ~ 0. The small error 
found in the calculation can be attributed to truncation errors in the finite difference 
approximations. These results match with analytical calculations. 
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CONVERGENCE GRAPH 




LEGEND 

□ = CONV FACTOR 



Figure 5. Convergence Graph 
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IV RESULTS AND CONCLUSIONS 



A. INITIAL CONDITIONS 



To properly apply the computer program generated, inputted data, or conditions, 
must be within the tolerances of the main equations used. This is demonstrated by the 
use of .01 meters as the thickness of the thin plate and a frequency range from 500 Hz 
to 4000 Hz. Other inputted data, including the densities of steel and water were taken 
directly from Kinsler, Frey, Coppens and Sanders [Ref. l:pp. 461-463]. 

The choice of the Ax and Ay terms is based on the need to have a significant 
number of data points in the waveguide domain that will give an accurate picture of what 
is happening. This is also true for the time step chosen, we need to make sure that the 
At was small enough to keep the problem from blowing up. A Courant-Lewy-Friedrichs 
stability condition must be satisfied for the domain which limits the size of At. 
Additionally as was mentioned before, the computer program was run for the same 
amount of time for each trial, this time period is long enough to ensure that a steady state 
solution was being achieved. 
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B. PROGRAM OUTPUT DATA 



The data obtained from the computer program demonstrates the following major 
points. When the frequency is increased you get more propagating modes from each 
particular plate length. When L is changed by a specific amount you do see a difference 
in the magnitude and shape of the deformation graphs. 

Our first table shows the calculated magnitudes of the propagating modes of the 
pressure wave at an assigned frequency of 750 Hz at different beam spacings. 

Table I. MAGNITUDE OF RADIATING MODES. 

For a frequency of 750 Hz 

with L = 0.8 fvo = .467808127 

with L= 1.0 cvo = -484098911 
a, = .307484090 

with L = 1.2 ao = .232444227 
a, = .367214561 



In addition to the data in Table I, we have in Figures 6 through 8 a graphical 
representation of the scalded pressure amplitude. The pressure was measured at the 
artificial boundary , yg at a frequency of 750 Hz. The horizontal axis is of a length of 
2L for each of the graphs starting at x = - 1 and going to x = 1 . The numbering on the 
horizontal axis represents each of the 41 pressure nodes used in the finite difference 
code. 
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SCALED PRESSURE AMPLITUDE 



L = 0.8 FREQUENCY = 750 HZ 




LEGEND 

Q = FLUID PRESSURE 



Figure 6. Sample Pressure Graph 



SCALED PRESSURE AMPUTUDE 



L = 1.0 FREQUENCY = 750 HZ 




LEGEND 

a = FLUID PRESSURE 



Figure 7. Sample Pressure Graph 



SCALED PRESSURE AMPLmJDE 



L = 1.2 FREQUENCY = 750 HZ 




LEGEND 

□ = FLUID PRESSURE 



Figure 8. Sample Pressure Graph 
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Table II is a list of the magnitudes of the a„’s for program runs with a frequency 



of 3500 Hz. 



Table II. MAGNITUDE OF RADIATING MODES. 



For a frequency of 3500 Hz 



with L = 0.8 


“o = 


.451181054 




a, = 


.323620617 




«2 = 


.088650584 


with L = 1.0 


“o = 


.365310729 




a-| = 


.260265648 




Ol2 = 


.062102031 




“3 = 


.060141488 




«4 = 


.076183497 


with L = 1.2 


CVo = 


.252027571 




a, = 


.165709317 




ao = 


.047427550 




«3 = 


.046469088 






.049672558 




05 = 


.043821129 




«6 = 


.067380726 



From the data above we obser\'e that at a frequency of 3500 Hz there is a 
substantial change in the magnitude of the propagating modes as L is increased from 0.8 
to 1.2. It can be seen from other runs that the higher the frequency the larger the 
differences in the magnitudes of the propagating modes. 

Figures 8 through 10 show the scaled pressure amplitude data obtained from the 
program runs with an assigned frequency of 3500 Hz. It is observed that changing of 
the length L on which the Ax is scaled has an effect on not only the amplitude of the 
pressure wave but also their shape. 
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SCALED PRESSURE AMPLITUDE 



L = 0.8 FREQUENCY = 3500 HZ 




LEGEND 

o = FLUID PRESSURE 



Figure 9. Sample Pressure Graph 
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L = 1.0 FREQUENCY = 3500 HZ 




LEGEND 

a = FLUID PRESSURE 



Figure 10. Sample Pressure Graph 



SCALED PRESSURE AMPLITUDE 



L = 1.2 FREQUENCY = 3500 HZ 




LEGEND 

D = FLUID PRESSURE 



Figure 11. Sample Pressure Graph 
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Table III is a list smaple energy calculations taken at three frequencies inputted into 
the computer program. The data presented demonstrates that for theses frequencies and 
lengths run in the computer program that the fundamental mode contains the largest 
amount of energy. Also evident is that the energy contained in the fundamental mode 
starts becomes more constant as the frequency is increased for different lengths of L. 
This data is in agreement with the theory of sound propagation. 
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Table HI ENERGY CALCULATIONS 



Percent of Energy in 


Percent of Energy in 


ProDaeatine Mode 


Prooaeatine Mode 


For a frequency of 750 Hz 


For a frequency of 3500 Hz 


with L = 0.8 Eo = 100% 


with L = 0.8 Eo = 66 % 
E, = 32% 


with L = 1.0 Eo = 99% 


E 2 = 2 % 


E, = 1 % 


with L = 1.2 Eo = 77% 


with L = 1.0 Eo = 65% 


E, = 23% 


El = 32% 

E 2 = 2 % 


For a frequency of 1750 Hz 


E 3 = 1 % 

E 4 = * 


with L = 0.8 Eo = 77% 


E, = 23% 


with L = 1.2 Eo = 64% 
E, = 27% 


with L = 1.0 Eo = 68 % 


E 2 = 2 % 


E, = 31% 


E 3 = 2 % 


Ei= 1% 


E 4 = 2 % 

E 5 = 1 % 


with L = 1.2 Eo = 64% 


E. = 2 % 


E, = 32% 
Ej = 3% 

E 3 = 1 % 

* denotes much less then 1 % 
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Table IV shows the magnitude of the primary propagating mode for each of the 



different frequencies and lengths inputted into the computer program. 
Table TV MAGNITUDES OF PRIMARY MODES 



Magnitude of ao’s 



L = 0.8 L = 1.0 



Frequency 
750 Hz 

1000 Hz 

1250 Hz 

1500 Hz 

1750 Hz 

2000 Hz 

2250 Hz 

2500 Hz 

2750 Hz 

3000 Hz 

3250 Hz 

3500 Hz 

3750 Hz 



.467808127 

.6182562721 

.553467214 

.498957574 

.404151559 

.486700475 

.394520938 

.385663450 

.385663450 

.420025572 

.369634628 

.451181054 

.446322470 



.484098911 

.444052100 

.457091272 

.415586054 

.431289613 

.433078885 

.400619507 

.392564356 

.392564354 

.353338659 

.226448834 

.365310729 

.330691218 



L = 1,2 
.367214561 
.262355387 
.343211770 
.320898831 
.338643909 
.316784620 
.252146780 
.267684340 
.267684340 
.234915197 
.194508076 
.252027571 
.226533830 
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The data obtained from the computer program and represented in Table III shows 
that there is always a difference in the magnitude of the fundamental propagating mode 
for changes in L at all frequencies comprised within the computer runs. The difference 
in the magnitude of the primary propagating mode between the runs with a length "L" 
of 1.2 and the runs with a length L of 1 are greater then the difference between the runs 
with lengths of 0.8 and 1.0. In the frequency range covered when L is increased from 
1.0 to 1.2 there is always a drop in the magnitude of the fundamental propagating mode. 
This is not always true when L is decreased to 0.8 from 1.0. Although the magnitude 
of these increases are smaller then the magnitude of the decreases this demonstrates that 
by increasing the spacing of the reinforcing beams will not always result in a decrease 
in the amplitude of the fundamental propagating mode in this frequency range. From this 
result we can conclude that additional factors play an important role in determining the 
magnitude of the propagating modes. 

Because the model is allowed to run until a steady state solution is reached there 
is no appreciable difference in the absolute pressure measured at the artificial boundary 
for any of the frequencies observed. However, it is evident that as the length of the plate 
"L” is increased there is a slight drop in the pressure measurement. This is true for all 
frequencies. Figure 12 is a sample plate deformation graph for L = 1.0 and a frequency 
of 1000 Hz. 

In conclusion, the computer program produces results as expected for the different 
frequencies of the plane wave and lengths of the thin plate. The superpositioning of the 
reflected plane wave has showed that as L is increased there is a larger number of 
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propagating modes accompanied by a change in the amplitude of each mode. Additional 
results could be explored for different angles of deflection as well as different plate 
thicknesses. This will require modification of the computer program as it stands to 
gather data for additional angles and thicknesses. 
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SCALED LENGTH OF 1.0 




LEGEND 

□ = DISP AT 1000 HZ 



Figure 12. Deformation Graph 
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APPENDIX: COMPUTER PROGRAM 



FILE. 280CT FORTRAN A1 



c XXXX )OOnOOO(X X X )000( X X X X X XX X XX X X X XX XXXX X K )00( X X X X X X X X)000( X X X xxx 

c 

. THIS PROGRAM SIMULATES THE DISTRUBANCE IN A FLUID 
AFTER A SOUND WAVE HAS HIT A PLATE. 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



P = PRESSURE (EXCESS) 

X = HORIZONTIAL POSITION 
W = DEFORMATION 
SPD = WAVE SPEED 
FREQ = FREQUENCY 



D = FLEXURAL RIGIDITY 
Y = VERTICAL POSITION 
KF = FLUID WAVE NUMBER 
KP = PLATE WAVE NUMBER 
OMEGA = RATIO OF (KF/KP)XX2 



COUNTERS. 



N 



THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER (IX = 20, lY = 80) 

PARAMETER (H = .025, DELTAT = .0325 ) 

PARAMETER (NDIM = 2XIX+1) 

THE FOLLOWING ARE COMPLEX ARRAYS 

0. lY, 5), W( -IX. IX, 3) 



P( -IX: IX, 

LB, FACR 

A( -IX. IX, -IX. 
AA( NDIM, NDIM) 
BB( NDIM, NDIM) 
C( -IX. IX) 

NDIM ) 



IX) 



X, PI, 



COMPLEXX16 
COMPLEX 
COMPLEX 
COMPLEX 
COMPLEX 
COMPLEX 

COMPLEX Z( 

C0MPLEXX16 FI 
COMMON/SS/P, W 

COMMON/TT/KP,KF,INFA, THETAI,L, 

1 EISP, OMEGA, T 

COMMON/CC/LB, FACR 
COMMON/EE/WGTWO, WGTHE, WGTOT, CON 
COMMON/PT/PGTWO, PGTHE, PGTOT 
COMMON/HH/MINN, MAXX 
COMMON/II/AA, BB 
COMMON/KK/IPVT 
COMMON/MM/NN 
COMMON/OO/C 
COMMON/FF/A 
COMMON/ UP/ A L 
COMMON/RN/RUNS 
COMMON/FEL/U,V,TP 
COMMON/SN/FREQ, SPD 
COMMON/GF/WGGRF, PGGRF 

INTEGER RUNS, IPVT(NDIM), SAM 
PARAMETER (SAM = 1000) 

REAL KF, KP, INFA, EISP, OMEGA, 

REALK8 WGTWO, WGTHE, WGTOT, CON 
REALK8 PGTWO, PGTHE, PGTOT 

REAL WGGRF( SAM ), PGGRF( SAM ) 

REAL THETAI, FREQ, SPD 

REAL X, T 

PI = 3.1A1592654 
CON = .015D0 



L,PI 
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no onoor 



WGTHO = O.ODO 
MGTHE = O.ODO 
PGTWO = O.ODO 
PGTHE = O.ODO 



CALL XUFLOW (0) 
CALL INIT 
RUNS = 0 
20 CALL FLPL 
CALL DOMAIN 

RUNS = RUNS + 1 



CALL ARTF 

IF (RUNS ,GT. 1999) THEN 

GOTO 60 

ELSE 

GOTO 61 

ENDIF 

60 CALL FINN 

61 T = T + DELTAT 
CALL SHIF 

IF (RUNS .EQ. 3000) THEN 

GOTO 80 

ELSE 

GOTO 20 

ENDIF 

80 CALL ALPHA 
CALL CLOSE 

WRITE(2,767)THETAI,KF 

767 FORMATdX, ’THETAI = ',F12.6,3X,' KF 
WRITE(2,768)OMEGA,EISP 

768 FORMATdX, ' OMEGA = • , F12 . 6 , 3X, ' EISP 
WRITE(2,196)DELTAT, KP 

196 FORMATdX, 'DELTAT = ',F12.6,3X,' KP 
HRITE(2,A87)H,L 

A87 FORMATdX,' H = ',F12.6,3X, ' L 
WRITE(2,735)FREQ,SPD 
735 FORMATdX, 'FREQ = ' , F12 . 6 , 3X, ' SOUND 
WRITE(2,629)PI 

629 FORMATdX,' PI = ',F12.7) 
WRITE(2,919) 

919 FORMATdX, 'FACR = ') 

MRITE(2,K) FACR 



= ',F12.6) 

= ',F12.7) 

= ',F12.6) 

= ',F6.3) 

SPEED = ',F12.6) 
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oooooo o ooooo 



WRITE(2,222) 

222 FORMATCIX, 'FINAL DATA* ,//, IX, ' X 

1 ' W ’ , • P ’ ) 

DO 333 I = -IX, IX, 1 
XX = I*H 

333 WRITE(2,666)XX,CDABS(W(I,2)),CDABS(P(I,IY,3)) 

666 FORMAT(1X,5F10.5) 

HRITE(2,223)RUNS 

223 > FORMATCIX, 'RUNS = ',17) 

WRITE(2,678) 

678 FORMATCIX, 'BELOW IS THE FINAL WGTOT, AND PGTOT') 
WRITEC2,x)WGT0T 
WRITEC2,K)PGT0T 
STOP 
END 

SUBROUTINE INIT 

THIS SETS UP THE INITIAL CONDITIONS OF THE 
PROBLEM INCLUDING THE DOMAIN AND THE 
ARTIFICIAL BOUNDARY 



C 



THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER CIX = 20, lY = 80) 

PARAMETER CH = .025, DELTAT = .0325) 

PARAMETER CNDIM = 2XIX+1) 

THE FOLLOWING ARE COMPLEX ARRAYS 



COMPLEXX16 

COMPLEX 

COMPLEX 

COMPLEX 

COMPLEX 



0. lY, 3), WC -IXi IX, 3) 



PC -IX! IX, 

LB, FACR 

AC -IX!lX, -IX! IX) 

AAC NDIM , NDIM 
ZC NDIM ) 

COMMON/SS/P, W 

COMMON/TT/KP,KF,INFA, THETAI,L, X, PI, 



), BBC NDIM 



NDIM 



I) 



FILE! 280CT FORTRAN A1 



C 



1 EISP, OMEGA, T 

COMMON/CC/LB, FACR 
COMMON/FF/A 

COMMON/GG/GAMMA, BETTA 
COMMON/ HH/MINN, MAXX 
COMMON/II/AA, BB 
COMMON/KK/IPVT 
COMMON/FEL/U,V,TP 
COMMON/SN/FREQ, SPD 



REAL THETAI, INFA, KF, EISP, OMEGA, 
REAL ART, K, XI, XIL, RCOND 
REAL GAMMAC-AO! 40), BETTA C-40!40) 

REAL SPD, FREQ 

REAL X, T 

INTEGER I, J, N, LJ, NN, U, V, TP 
INTEGER NMODMIN, NMODMAX, MINN, MAXX 
INTEGER IPVTCNDIM) 



TIN, 

PI 



L, KP 
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ooo o o oo 



TP = <(1 

LB = CMPLX(0.0, 1.0) 

FACR = CMPLXd.O, 0.0) 

SPD = 1500. ODO 
FREQ = 3750. ODO 
N = 0 
1 = 0 
J = 0 
NH = 1 
L = 1.2D0 

THETAI = PI/2. ODO 

T = DELTAT 
NMODMIN = 0 
NMODMAX = 0 

KF = (2D0XPIXFREQXD/SPD 

. GT = ((195000000000. 0*( .05XX3))/12(l-( .28XX2)) 
KP = ((((FREQX2»PI)XX2)X7700XH)/FR)XX.25 
KP = A0.75349751A7D0 
OMEGA = (KF/KP)XX2 

INFA = (DELTATXX2)/((KFXX2)X(LXX2)X(HXX2)) 

DO 12 I = -IX,IX,1 
DO 11 J = 0,IY»1 

P(I,J,1) = CMPLX(0.0, 0.0) 

P(I,J,2) = CMPLX(0.0, 0.0) 

P(I,J,3) = CMPLX(0.0, 0.0) 

11 CONTINUE 

12 CONTINUE 

DO 21 I = -IX, IX, 1 

H(I,1) = CMPLX(0.0, 0.0) 

H(I,2) = CMPLX(0.0, 0.0) 

H(I,3) = CMPLX(0.0, 0.0) 

21 CONTINUE 

THIS SETS UP THE ARTIFICAL BOUNDRAY 
K = KFXL 

TIN = -KXCOS(THETAI) 

DO 15 N = -40, 40, 1 

GAMMA(N) = TIN + PIXN 

IF ( ABS(GAHMA(N)) .LE. K) THEN 
NMODMIN = MIN(N, NMODMIN) 

NMODMAX = MAX(N, NMODMAX) 

ENDIF 

15 CONTINUE 

MINN = NMODMIN 
MAXX = NMODMAX 



DO 18 N = MINN, MAXX, 1 

BETTA(N) = SQRT((KXX2) - (GAMMA(N)XX2)) 
18 CONTINUE 

DO 78 I = -IX, IX, 1 
DO 79 LJ = -IX, IX, 1 

A(I,LJ) = CMPLX(0. 0,0.0) 

79 CONTINUE 
78 CONTINUE 

ART = (2X(HX»2))/(4XDELTAT) 

DO 300 I = -IX, IX, 1 

DO 301 LJ = -IX, IX, 1 
XI = IXH 
XIL = LJXH 

IF (ABS(LJ) .EQ. IX) THEN 
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DEL = .5 
ELSE 
DEL = 1 
ENDIF 

DO 302 N = MINN, MAXX, 1 
A(I,LJ) = A(I,LJ) + ARTKDEL 

1 KBETTA(N)XCEXP(LBKGAMMA(N)X(XI-XID) 

302 CONTINUE 
301 CONTINUE 
300 CONTINUE 
C 

DO A20 U = 1, NDIM, 1 
> DO A21 V = 1, NDIM, 1 
I = -IX+U-1 
LJ = -IX+V-1 

AA(U,V) = INFAX(A(I,LJ)) 

BB(U,V) = INFA*(A(I,LJ)) 

IF (I .EQ. LJ) THEN 

AA(U,V) = AA(U,V) + 1.0 
BB(U,V) = BB(U,V) - 1.0 

ENDIF 

A21 CONTINUE 
A20 CONTINUE 

CALL CGECO (AA, NDIM, NDIM, IPVT, RCOND,Z) 

RETURN 

END 

C 

SUBROUTINE FLPL 

THIS SUBROUTINE DEFINES THE FLUID-PLATE INTERFACE 
THIS IS WHERE Y=0 ALONG THE TOP OF THE PLATE. 

THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER (IX = 20, lY = 80) 

PARAMETER (H = .025, DELTAT = .0325) 

PARAMETER (NDIM = 2XIX+1) 

THE FOLLOWING ARE COMPLEX ARRAYS 
COMPLEXX16 P( -IX: IX, 0: lY, 3), W( -IX: IX, 3) 
COMPLEX LB, FACR, PTEMP 

C0MPLEXX16 FI 
COMMON/SS/P, W 

COMMON/TT/KP,KF,INFA, THETAI,L, X, PI, 

1 EISP, OMEGA, T 

COMMGN/CC/LB, FACR 
COMMON/RN/RUNS 

REAL KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
REAL X, T 

REALXA SQ, SW, BWl , BW2, BW3, BWA 
INTEGER I, J, RUNS 



EISP = 0.1310186A566D0 
J = 0 

SQ = (DELTAT/(H*KF))XX2 
SW = (DELTAT/(HxkP)kx2)kx 2 
BWl = 2-6?(SW 
BH2 = AXSW 

BW3 = EISPKDELTATXK2XKF/0MEGA 
BWA = -SW 

A = FL0AT(IX)/FL0AT(2XIX) 
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THIS LOOP COMPUTES THE INTERIOR NODES 

DO 10 I = -IX + 2, IX -2, 1 
X = IXH 

W(I,3) = BW1XW(I,2) -W(I,1) + BW2x(W(I+l,2) + Wa-l,2)) 

1 +BW4*(W(I+2,2)+W(I-2,2)) 

2 +BW3*(P(I,0,2)+ F1(X,T)) 

PTEMP = -(2*H/DELTATxx2)x(W(I,3)-2*W(I,2)+W(I,l)) +P(I,1,2) 

P(I,J,3) = INFA*(P(I+1, J,2)+P(I-1, J,2)-A*P(I, J,2)+P(I, J+l,2)+ 

1 PTEMP )+2xP(I,J,2) -P(I,J,1) 

10 CONTINUE 

THESE LOOPS COMPUTE THE NODES AT -19 AND 19 

I = -IX+1 
X “ ( ”IX+ 1 ) 

W(I,3) = BW1XW(I,2)-N(I,1) +BW2X(W(I+1,2)+W(I-1,2)) 

1 +BWA*(W(I+2,2)+W(I,2)) 

2 +BW3*(P(I,0,2)+ F1(X,T)) 

PTEMP = -(2*H/DELTAT)f)(2)x(W(I,3)-2xW(I,2)+W(I,l)) +P(I,1,2) 

P(I,J,3) = INFA*(P(I+1, J,2)+P(I-1,J,2)-AXP(I, J,2)+P(I, J+l,2)+ 

1 PTEMP )+2*P(I,J,2) -P(I,J,1) 



I = IX -1 
X = (IX-l)XH 

N(I,3) = BHl?(Na,2) - H(I,1) + BW2X( W( I+l , 2)+H( I-l , 2) ) 

1 +BUAx(W(I,2)+W(I-2,2)) 

2 +BW3X(P(I,0,2)+ F1(X,T)) 

PTEMP = -(2XH/DELTAT*X2)X(W(I,3)-2XW(I,2)+W(I,1)) +P(I,1,2) 
P(I,J,3) = INFAX(P(I+1,J,2)+P(I-1,J,2)-AXP(I,J,2)+P(I, J+l,2)+ 

1 PTEMP )+2*P(I,J,2) -P(I,J,1) 

THESE ARE THE END POINTS 

FACR=CEXP(2xLBxLxKFxC0S(THETAI) ) 

I = -IX 
J = 0 

W(I,3) = CMPLXIO. 0,0.0) 

PTEMP = P(I,1,2) 

P(I,J,3) = INFAX(P(I+1,J,2)+FACRXP(-I-1, J,2)-AXP(I, J,2) 

1 +P(I, J+1,2)+PTEMP )+2xPCI,J,2) -P(I,J,1) 

FACR=CEXP(-2XLB*LXKFxC0S(THETAI)) 

I = IX 
J = 0 

H(I,3) = CMPLXIO. 0,0.0) 

PTEMP = P(I,1,2) 

P(I,J,3) = INFAX(FACRXP(-I+1, J,2)+P(I-1, J,2)-A*P(I,J,2) 

1 +P(I, J+1,2)+PTEMP )+2XP(I,J,2) -P(I,J,1) 

RETURN 

END 
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FORTRAN A1 



SUBROUTINE DOMAIN 

THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER (IX = 20, lY = 80) 

PARAMETER (H = .025, DELTAT = .0325) 

PARAMETER (NDIM = 2XIX+I) 

THE FOLLOWING ARE COMPLEX ARRAYS 
C0MPLEXX16 P( -IXt IX, 0. lY, 3), W( -IXt IX, 3) 
COMPLEX LB, FACR 

COMMON/SS/P, W 

COMMON/TT/KP, KF, INFA, THETAI, L, X, PI, 

I EISP, OMEGA ,T 

COMMON/CC/LB, FACR 

REAL KP, KF, INFA, THETAI, L, PI, EISP, OMEGA 
REAL X, T 
INTEGER I,J 

DO 22 I = I-IX, IX-1,I 
DO 21 J = 1, IY-1, 1 

P(I,J,3) = INFAX(P(I+1, J,2)+P(I-1, J,2)- 

1 AXP(I,J,2)+P(I, J+1,2)+P(I, J-l,2))+ 

2 2XP(I, J,2)-P(I,J,1) 

21 CONTINUE 

22 CONTINUE 

LB = CMPLX (0.0, 1,0) 

LA = LKKFxCOS(THETAI) 

FACR = CEXP((2)X(LBXLA)) 

I = -IX 

DO A1 J = 1, IY-1 , 1 

P(I,J,3) = INFAX(P(I+1, J,2)+ FACRKP(-I-1, J,2) 

1 -AXP(I,J,2)+P(I, J+1,2)+P(I,J-1,2))+ 

2 2XP(I,J,2)-P(I, J,l) 

A1 CONTINUE 

LB = CMPLX (0.0, 1.0) 

RA = L)(KFXCOS(THETAI) 

FACR = CEXP((-2)X(LBXRA)) 

I = IX 

DO 31 J = 1, IY-1,1 

P(I,J,3) = INFAK(FACRXP(-I+1, J,2)+ P(I-1,J,2) 

1 -A«P(I, J,2)+P(I, J+1,2)+P(I,J-1,2))+ 

2 2XP(I, J,2)-P(I, J,l) 

31 CONTINUE 

RETURN 

END 
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SUBROUTINE ARTF 

THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER (IX = 20, lY = 80) 

‘PARAMETER (H = .025, DELTAT = .0325) 

PARAMETER (NDIM = 2XIX+1) 

THE FOLLOWING ARE COMPLEX ARRAYS 



PC -IX« IX, Ot lY, 3), M( -IXi IX, 3) 

LB, FACR 

AC *"IX* IX “IX* IX) 

AA( NDIM ', NDIM ), BBC NDIM , NDIM 
BCNDIM) 

CC -IX: IX) 

) 



C0MPLEXX16 
COMPLEX 
COMPLEX 
COMPLEX 
COMPLEX 
COMPLEX 

COMPLEX Z( NDIM 

COMMON/SS/P, W 

COMMON/TT/KP,KF,INFA, THETAI,L,X, PI, 
1 EISP, OMEGA ,T 

COMMON/ CC/LB, FACR 
COMMON/ FF/ A 
COMMON/II/AA, BB 



SUBROUTINE SHIF 

THIS MOVES THE PROGRAM TO THE NEXT TIME LEVEL 



THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER (IX = 20, lY = 80) 

PARAMETER (H = .025, DELTAT = .0325) 

PARAMETER (NDIM = 2XIX+1) 

C THE FOLLOWING ARE COMPLEX ARRAYS 

COMPLEX«16 P( -IX: IX, 0: lY, 3), W( -IX: IX, 3) 
COMPLEX LB, FACR 

COMMON/SS/P, W 

COMMON/TT/KP,KF,INFA, THETAI, L, X, PI, 

I EISP, OMEGA, T 

COMMON/CC/LB, FACR 
C 

REAL KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
REAL X, T 
INTEGER I, J 
C 
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C0MM0N/FEL/U,V,TP 

tOMMON/KK/IPVT 

C0MM0^i/00/C 

REAL KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
REAL X, T 

INTEGER IPVT(NDIM), JOB ,U,V, TP 
TP = 

JOB = 0 

LB = CMPLX (0.0, I .0) 

LA = LkKFxcOS(THETAI) 

RA = LXKFXCOS(THETAI) 

DO A2 I = -IX+1, IX-1, I 

C(I) = INFAX(P(I+1,IY,2)+P(I-1,IY,2)+2XP(I,IY-1,2)) 

1 + (2 - (AXINFA))XP(I,IY,2) 

A2 CONTINUE 
I = -IX 

FACR = CEXP((2)X(LBXLA)) 

C(I) = INFAX(P(I+1,IY,2)+FACRXP(-I-1,IY,2)+2XP(I,IY-1,2)) 
1 + (2 - (AXINFA))XP(I,IY,2) 

I = IX 

FACR = CEXP((-2)X(LBXRA)) 

C(I) = INFAX(FACRXP(I-I,IY,2)+P(-I+1,IY,2)+2XP(I,IY-1,2)) 
I + (2 - (AXINFA))XP(I, IY,2) 

DO 10 I = 1, NDIM, 1 

BCD = INFAX(CMPLX(0.0, 0.0)) 

10 CONTINUE 



DO 20 I = 1,NDIM,1 

DO 30 LJ = 1, NDIM, 1 

B(I) = B(I) + BB(I,LJ)XP(-IX+LJ-1,IY,1) 
30 CONTINUE 
20 CONTINUE 

DO 150 I = 1, NDIM, 1 

B(I) = B(I) + C(-IX+I-1) 

150 CONTINUE 

CALL CGESLCAA, NDIM, NDIM, IPVT,B, JOB) 

DO A6 I = -IX, IX, 1 

P(I,IY,3) = B(IX+1+I) 

A6 CONTINUE 
RETURN 
END 
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DO 19 I = - 


IX, 


IX, 


1 




DO 22 J = 


0, 


lY, 


1 




P(I, J,l) 


= P(I,s 


J,2) 


22 


CONTINUE 








19 


CONTINUE 










DO 20 I = - 


IX, 


IX, 


1 




DO 23 J = 


0, 


lY, 


1 




P(I,J,2) 


= 1 


P(I,. 


J,3) 


23 


CONTINUE 








20 


CONTINUE 










DO 29 I = - 


IX, 


IX, 


1 




W(I,1) = 


MC 


1,2) 




29 


CONTINUE 










DO 90 I = - 


IX, 


IX, 


1 




W(I,2) = 


HC 


1,3) 




90 


CONTINUE 









RETURN 

END 



FUNCTION FKXP, TIME) 

PARAMETER (IX = 20, lY = 80) 

PARAMETER (H = .025, DELTAT = .0325) 

C0MPLEXXI6 P( -IX! IX, 0i lY, 5), W( -IXi IX, 3) 
C0MPLEXXI6 FI, SI 
COMMON/SS/P, W 
COMMON/RN/RUNS 

COMMON/TT/KP,KF,INFA, THETAI,L, X, PI, 

I EISP, OMEGA, T 

REAL KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
REAL X, T 
REALMS Q9,Q5,S 



IF(RUNS .EQ. 1) THEN 

S = 0.0 

ELSE 

S = TIME+XP)(KFXCOS(THETAI) 
ENDIF 
Q9 = 0 

Q5 = AMINKTIME, 18.0) 
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Q8 = DEXP(2XQ5-6 .0) 
Q7 = Q8/(Q8+,1/Q8) 
SI = CMPLX(Q9,-S) 

FI = Q7KCDEXP(S1) 

RETURN 

END 



, SUBROUTINE ALPHA 



THIS SETS UP THE INITIAL CONDITIONS OF THE 
PROBLEM INCLUDING THE DOMAIN AND THE 
ARTIFICIAL BOUNDARY 



THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER (IX = 20, lY = 80) 

PARAMETER (H = .025, DELTAT = .0325) 

PARAMETER (NDIM = 2XIX+1) 

THE FOLLOWING ARE COMPLEX ARRAYS 
C0MPLEXX16 P( -IX. IX, 0. lY, 3), W( -IX. IX, 3) 
COMPLEX LB, FACR 

COMPLEX A( -IX. IX, -IX. IX) 

COMPLEX AA( NDIM , NDIM ), BB( NDIM , NDIM 

COMPLEX AL( -AO. AO) 



COMPLEX LAND 

COMMON/SS/P, W 

COMMON/TT/KP,KF,INFA, THETAI,L, X, PI, 
1 EISP, OMEGA, T 

COMMON/CC/LB, FACR 
COMMON/FF/A 

COMMON/GG/GAMMA, BETTA 
COMMON/HH/MINN, MAXX 
COMMON/II/AA, BB 
COMMON/KK/IPVT 
COMMON/RN/RUNS 
COMMON/FEL/U,V,TP 
COMMON/SN/FREQ, SPD 



REAL THETAI, INFA, KF, EISP, OMEGA, L, KP 
REAL XI, X, T 

REAL GAMMA(-A0. AO), BETTA (-A0.A0), PI* 

REAL SPD, FREQ 

INTEGER I, U, V, TP, MINN, MAXX, IPVT(NDIM), RUNS, Y 



Y = IY/2 
WRITE(2,75) 

75 FORMATdX,' ') 
HEN = RUNSKDELTAT 



DO 10 M = MINN, MAXX,1 

AL(M) = CMPLX(0.0, 0.0) 
DO 30 I = -IX, IX, 1 
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XI = IXH 

IF( ABS(I) .EQ. IX) THEN 
DEL = .5 
ELSE 

DEL = 1.0 
ENDIF 

LAND = CEXP(-LBXYXHXBETTA(M) + (LBxHEN)) 

AL(M) = AL(M) + .5XLANDXHXDELXCEXP(-LBXGAMHA(M)xXI)xP(I,Y,5) 
30 CONTINUE 
10 CONTINUE 
WRITE(2,32) 

32 FORMATUX, 'ALPHA •’N"S FOR RADIATING MODES') 

DO 111 M = MINN, MAXX 
HRITE(2,)()CABS(AL(M)) 

111 CONTINUE 
RETURN 
END 

SUBROUTINE CLOSE 

THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 
PARAMETER (IX = 20, lY = 80) 

PARAMETER (H = .025, DELTAT = .0325) 

PARAMETER (XYX= AO) 

PARAMETER (TNT= 1000) 

PARAMETER (NDIM = 2XIX+1) 

C THE FOLLOWING ARE COMPLEX ARRAYS 

C0MPLEXX16 P( -IX: IX, 0: lY, 3), W( -IXi IX, 3) 

COMPLEX LB, FACR 

C THE VALUES IN THE COMMON "BLOCKS" ARE VARIABLES 

COMMON/SS/P, W 

COMMON/TT/KP,KF, INFA, THETAI,L, X, PI, 

1 EISP, OMEGA ,T 

COMMON/CC/LB, FACR 
COMMON/HH/MINN, MAXX 
COMMON/GF/WGGRF, PGGRF 
C 

REAL KP, KF, INFA, THETAI, L, PI, EISP, OMEGA 



REAL X, T 

REAL XX( 0. Al), WPLOT( -IX: IX), ZZ(O-.IOOO) 

REAL PPLOSC -IX: IX), SMALL(IOOO) 

INTEGER II, RUNS, SAM, IR 
PARAMETER (SAM = 1000) 

REAL WGGRF( SAM ), PGGRF( SAM ) 

C 

DO 10 II = 0, XYX, 1 
XX(II) = II 

10 CONTINUE 
C 

DO 11 IR = 0, TNT, 1 
ZZ(IR) = IR 

11 CONTINUE 
C 

DO 30 I = 1,1000,1 

SMALL! I) = WGGRF(I) 

30 CONTINUE 

DO 35 I = -IX, IX, 1 

PPLOS(I) = CDABS(P(I,IY,3)) 

35 CONTINUE 

DO A5 I = -IX, IX, 1 

WPLOT(I) = CDABS(M(I,3)) 

A5 CONTINUE 
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CALL COMPRS 

CALL PLOTDCZZ, SMALL, 1000, .TRUE. , ' LINLIN* , 'CONV FACTOR^’, 

1 'CONVERGENCE GRAPHS', 'RUN NUMBER + 2000$', 

2 'CONVERGENCE VALUED') 

CALL PLOTD(XX,PPLOS,Al, .TRUE. ,' LINLIN' ,' FLUID PRESSURES ' , 

1 'PI/2 LENGTH = 1.2 FREQ= 3750$ ', 'HORIZONTAL POSITION? ' , 

2 'SCALED PRESSURE?') 

CALL PLOTD(XX,WPLOT, Al, .TRUE. , ' L INLIN' , ' DISP AT 3750 HZ?', 

1 'SCALED LENGTH OF 1 .2? ', 'SCALED HORIZONTAL POSITION?', 

2 'SCALED DISPLACEMENT OF PLATE?') 

CALL DONEPL 

RETURN 
, END 

SUBROUTINE FINN 

THIS SUBROUTINE MAKES SURE WE HAVE REACHED THE STEADY STATE 
SOLUTION OF THE PROBLEM, IE ' EOUALIBRIUM' 

THE VALUES IN THE PARAMETER STATEMENTS ARE CONSTANTS 

PARAMETER (IX = 20, lY = 80) 

PARAMETER (H = .025, DELTAT = .0325) 

THE FOLLOWING ARE COMPLEX ARRAYS 
C0MPLEXX16 P( -IXi IX, 0. lY, 3), W( -IX. IX, 3) 

COMMON/SS/P, W 

COMMON/TT/KP,KF,INFA, THETAI, L, X, PI, 

1 EISP, OMEGA ,T 

COMMON/EE/WGTWO, WGTHE, WGTOT, CON 
COMMON/PT/PGTWO, PGTHE, PGTOT 
COMMON/RN/RUNS 
COMMON/GF/WGGRF, PGGRF 

REAL KP, L, OMEGA, EISP, KF, INFA, PI, THETAI 
REAL X, T 

REALX8 WGTWO, WGTHE, WGTOT, CON 
REALK8 PGTWO, PGTHE, PGTOT 
INTEGER I, RUNS, SAM 
PARAMETER (SAM = 1000) 

REAL WGGRF( SAM ), PGGRF( SAM ) 

DO 3A I = -IX+2, IX-2, 1 



WGTWO = W(I,2)KDC0NJG(W(I,2))+ WGTWO 

WGTHE = W(I,3)*DC0NJG(W(I,3))+ WGTHE 

WGTOT = (WGTHE - WGTWO) / WGTHE 

3A COflTINUE 



DO A4 I = -IX , IX , 1 

PGTWO = P(I,IY,2)XDC0NJG(P(I,IY,2))+ PGTWO 

PGTHE = P(I,IY,3)xdC0NJG(P(I,IY,3))+ PGTHE 

PGTOT = . (PGTHE - PGTWO) / PGTHE 

4A CONTINUE 



IF (RUNS.GT. 1999) THEN 
WGGRF(RUNS-1999) = WGTOT 
PGGRF(RUNS-1999) = PGTOT 
ENDIF 



WRITE(2,88) 

88 FORMATdX, 'WGGRF, PGGRF IN FINN') 
WRITE(2,X) WGGRF(RUNS-1999) 
WRITE(2,X) PGGRF(RUNS-1999) 

RETURN 

END 
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